{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoscilabs.inversion.LinearInversionDirect import LinearInversionDirectApp\n",
    "from ipywidgets import interact, FloatSlider, ToggleButtons, IntSlider, FloatText, IntText"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "app = LinearInversionDirectApp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Inversion App\n",
    "\n",
    "This app is based upon the inversion tutorial: \"INVERSION FOR APPLIED GEOPHYSICS\" by Oldenburg and Li (2005).\n",
    "\n",
    "Douglas W. Oldenburg and Yaoguo Li (2005) 5. Inversion for Applied Geophysics: A Tutorial. Near-Surface Geophysics: pp. 89-150.\n",
    "eISBN: 978-1-56080-171-9 \n",
    "print ISBN: 978-1-56080-130-6 \n",
    "https://doi.org/10.1190/1.9781560801719.ch5 \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Purpose\n",
    "\n",
    "By using a simple decaying and oscillating kernel function, which emulates the physics of electromagnetic (EM) survey, we understand basic concepts of inverting data. Three items that we are going to explore are:\n",
    "\n",
    "- Step1: Generate a sensitivity kernel (or matrix), $\\mathbf{G}$\n",
    "- Step2: Set a model then Simulate data ($\\mathbf{d} = \\mathbf{G}\\mathbf{m}$)\n",
    "- Step3: Invert the data, and explore inversion results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward problem\n",
    "\n",
    "\n",
    "Let $g_j(x)$ denote the kernel function for $j$th datum. With a given model $m(x)$, $j$th datum can be computed by solving following integral equation:\n",
    "\n",
    " $$ d_j = \\int_a^{b} g_j(x) m(x) dx $$\n",
    "\n",
    "where \n",
    "\n",
    "$$ g_j(x) = e^{jpx} cos (2 \\pi jqx) $$ \n",
    "\n",
    "By discretizing $g_j(x)$ we obtain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\mathbf{g}_j(\\mathbf{x}) = e^{jp\\mathbf{x}} cos (2 \\pi jq \\mathbf{x})$$\n",
    "\n",
    "where\n",
    "\n",
    "- $\\mathbf{g}_j$: $j$th row vector for the sensitivty matrix ($1 \\times M$)\n",
    "- $\\mathbf{x}$: model location ($1 \\times M$)\n",
    "- $j_k$: k=1, 2, ..., N\n",
    "- $p$: decaying constant (<0)\n",
    "- $q$: oscillating constant (>0)\n",
    "\n",
    "By stacking multiple rows of $\\mathbf{g}_j$, we obtain sensitivity matrix, $\\mathbf{G}$: \n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbf{G} = \n",
    "    \\begin{bmatrix}\n",
    "        \\mathbf{g}_1\\\\\n",
    "        \\vdots\\\\\n",
    "        \\mathbf{g}_{N}\n",
    "    \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "Here, the size of the matrix $\\mathbf{G}$ is $(N \\times M)$. \n",
    "Finally data, $\\mathbf{d}$, can be written as a linear equation:\n",
    "\n",
    "$$ \\mathbf{d} = \\mathbf{G}\\mathbf{m}$$\n",
    "\n",
    "where $\\mathbf{m}$ is an inversion model; this is a column vector ($M \\times 1$). \n",
    "\n",
    "In real measurments, there will be various noises source, and hence observation, $\\mathbf{d}^{obs}$, can be written as \n",
    "\n",
    "$$ \\mathbf{d}^{obs} = \\mathbf{G}\\mathbf{m} + \\mathbf{noise}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Step1: Generate a sensitivity kernel (or matrix), $\\mathbf{G}$\n",
    "\n",
    "By using the following app, we explore each row vector of the sensitivity matrix, $\\mathbf{g}_j$. Parameters of the apps are:\n",
    "\n",
    "- M: # of model parameters\n",
    "- N: # of data\n",
    "- p: decaying constant (<0)\n",
    "- q: oscillating constant (>0)\n",
    "- ymin: maximum limit for y-axis\n",
    "- ymax: minimum limit for y-axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5a9710f39d1f4bb69001494e0ecdcdec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(IntSlider(value=20, continuous_update=False, description='N', min=1), IntSlider(value=100, continuous_update=False, description='M', min=1), FloatSlider(value=-0.15, continuous_update=False, description='p', max=0.0, min=-1.0, step=0.05), FloatSlider(value=0.25, continuous_update=False, description='q', max=1.0, step=0.05), FloatText(value=1.0, description='j1'), FloatText(value=19.0, description='jn'), ToggleButtons(description='scale', index=1, options=('linear', 'log'), value='log'), Checkbox(value=False, description='fixed'), FloatText(value=-0.005, description='ymin'), FloatText(value=0.011, description='ymax'), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q1 = app.interact_plot_G()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step2: Generate a model and simulate data\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model $m$ is a function defined on the interval (-2,2). Here we generate a model that is the sum of a: (a) background $m_{ref}$, (b) box car $m_1$ and (c) Gaussian $m_2$. The box car is defined by\n",
    "- m$_{background}$ : amplitude of the background\n",
    "- m1 : amplitude\n",
    "- $m1_{center}$ : center\n",
    "- $m1_{width}$ : width\n",
    "the Gaussian is defined by \n",
    "- m2 : amplitude\n",
    "- $m2_{center}$ : center\n",
    "- $m2_{sigma}$ : width of Gaussian (as defined by a standard deviation)\n",
    "\n",
    "The $j$th datum is the inner product of the $jth$ kernel $g_j(x)$ and the model $m(x)$. In discrete form it can be written as the dot product of the vector $g_j$ and the model vector $m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $$ d_j = \\mathbf{g}_j \\mathbf{m} $$\n",
    "\n",
    "If there are $N$ data, these data can be written as a column vector, $\\mathbf{d}$:\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbf{d} = \\mathbf{G}\\mathbf{m} = \n",
    "    \\begin{bmatrix}\n",
    "        d_1\\\\\n",
    "        \\vdots\\\\\n",
    "        d_{N}\n",
    "    \\end{bmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding Noise\n",
    "\n",
    "Observational data are always contaminated with noise. Here we add Gaussian noise $N(0,\\epsilon)$ (zero mean and standard deviation $\\sigma$). Here we choose \n",
    "\n",
    "$$ \\epsilon = \\% |d| + \\text{floor} $$\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b84511884d2415c8e6d2f6bdcd6f7ea",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, continuous_update=False, description='m$_{background}$', max=2.0, min=-2.0, step=0.05), FloatSlider(value=1.0, continuous_update=False, description='m1', max=2.0, min=-2.0, step=0.05), FloatSlider(value=0.2, continuous_update=False, description='m1$_{center}$', max=2.0, min=-2.0, step=0.05), FloatSlider(value=0.2, continuous_update=False, description='m1$_{width}$', max=0.5, step=0.05), FloatSlider(value=2.0, continuous_update=False, description='m2', max=2.0, min=-2.0, step=0.05), FloatSlider(value=0.75, continuous_update=False, description='m2$_{center}$', max=2.0, min=-2.0, step=0.05), FloatSlider(value=0.07, continuous_update=False, description='m2$_{sigma}$', max=0.1, min=0.01, step=0.01), SelectMultiple(description='option', index=(1,), options=('kernel', 'model', 'data'), value=('model',)), Checkbox(value=True, description='add_noise'), FloatText(value=5.0, description='percentage'), FloatText(value=0.02, description='floor'), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q2 = app.interact_plot_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inverse Problem\n",
    "\n",
    "In the inverse problem we attempt to find the model $\\mathbf{m}$ that gave rise to the observational data $\\mathbf{d}^{obs}$. The inverse problem is formulated as an optimization problem: \n",
    "\n",
    "\n",
    "$$\\text{minimize} \\ \\ \\ \\phi(\\mathbf{m}) = \\phi_d(\\mathbf{m}) + \\beta \\phi_m(\\mathbf{m}) $$\n",
    "\n",
    "where \n",
    "\n",
    "- $\\phi_d$: data misfit\n",
    "- $\\phi_m$: model regularization\n",
    "- $\\beta$: trade-off (or Tikhonov) parameter  $0<\\beta<\\infty$\n",
    "\n",
    "Data misfit is defined as \n",
    "\n",
    "$$ \\phi_d = \\sum_{j=1}^{N}\\Big(\\frac{\\mathbf{g}_j\\mathbf{m}-d^{obs}_j}{\\epsilon_j}\\Big)^2$$\n",
    "\n",
    "where $\\epsilon_j$  is an estimate of the standard deviation of the $j$th datum.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model regularization term, $\\phi_m$, can be written as \n",
    "\n",
    "$$ \\phi_m(\\mathbf{m}) = \\alpha_s \\int (\\mathbf{m}-\\mathbf{m}_{ref}) dx + \\alpha_x \\int (\\frac{d \\mathbf{m}}{dx}) dx$$\n",
    "\n",
    "The first term is referred to as the \"smallness\" term. Minimizing this generates a model that is close to a reference model $m_{ref}$. The second term penalizes roughness of the model. It is generically referred to as a \"flattest\" or \"smoothness\" term.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step3: Invert the data, and explore inversion results\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the inverse problem we define parameters needed to evaluate the data misfit and the model regularization terms. We then deal with parameters associated with the inversion.\n",
    "\n",
    "### Parameters\n",
    "\n",
    "- `mode`: `Run` or `Explore`\n",
    "    - `Run`: Each click of the app, will run `n_beta` times of inversion\n",
    "    - `Explore`: Not running inversions, but explore result of the inversions\n",
    "\n",
    "- `mref`: reference model\n",
    "\n",
    "- `percent`: percentage of the uncertainty (%)\n",
    "\n",
    "- `floor`: floor of the uncertainty (%)\n",
    "\n",
    "- `chifact`: chi factor for stopping criteria (when $\\phi_d^{\\ast}=N \\rightarrow$ `chifact=1`)\n",
    "\n",
    "- `data`: `obs/pred` or `misfit`\n",
    "    - `obs/pred`: show observed and predicted data\n",
    "    - `misfit`: show normalized misfit\n",
    "\n",
    "- `beta_min`: minimum $\\beta$\n",
    "\n",
    "- `beta_max`: maximum $\\beta$\n",
    "\n",
    "- `n_beta`: the number of $\\beta$\n",
    "\n",
    "- `alpha_s`: $\\alpha_s$ for smallness\n",
    "\n",
    "- `alpha_x`: $\\alpha_x$ for smoothness\n",
    "\n",
    "- `option`: `misfit` or `tikhonov`\n",
    "    - `misfit`: show $\\phi_d$ and $\\phi_m$ as a function of $\\beta$\n",
    "    - `tikhonov`: show tikhonov curve\n",
    "    \n",
    "- `i_beta`: i-th $\\beta$ value\n",
    "\n",
    "- `scale`: `linear` or `log`\n",
    "    - `linear`: linear scale for plotting the third panel\n",
    "    - `log`: log scale for plotting the third panel     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e4799f7cd2724a8b80ca4bd69a0b30bd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(RadioButtons(description='mode', options=('Run', 'Explore'), value='Run'), FloatSlider(value=0.0, continuous_update=False, description='mref', max=2.0, min=-2.0, step=0.05), FloatText(value=5.0, description='percentage'), FloatText(value=0.02, description='floor'), FloatText(value=1.0, description='chifact'), ToggleButtons(description='data', options=('obs/pred', 'misfit'), value='obs/pred'), FloatText(value=0.001, description='beta_min'), FloatText(value=100000.0, description='beta_max'), IntText(value=81, description='n_beta'), FloatText(value=1.0, description='alpha_s'), FloatText(value=0.0, description='alpha_x'), ToggleButtons(description='option', index=1, options=('misfit', 'tikhonov'), value='tikhonov'), IntSlider(value=0, continuous_update=False, description='i_beta', max=80), ToggleButtons(description='scale', index=1, options=('linear', 'log'), value='log'), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q3 = app.interact_plot_inversion()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
